# -*- coding: utf-8 -*-
"""
Created on Sun Jun 30 13:48:38 2013

@author: Maxim
"""

import wind_tunnel_comparison as wt
import numpy as np
import matplotlib.pyplot as plt

def check_at_fc():
    pass

def mrc_translate(refOld,refNew,refLen,Cm,CL,CD,alpha):
    alpha = np.radians(alpha)
    CN = CL*np.cos(alpha) + CD*np.sin(alpha)
    CT = CD*np.cos(alpha) - CL*np.sin(alpha)
    CmNew = Cm + CN*(refNew[0]-refOld[0])/refLen - CT*(refNew[2]-refOld[2])/refLen
    return CmNew

def get_deriv(x,y):
    assert len(x)==len(y)
    deriv = np.zeros(len(x))
    deriv[0] = (y[1]-y[0])/(x[1]-x[0])
    deriv[-1] = (y[-1]-y[-2])/(x[-1]-x[-2])
    for i in range(len(x)-2):
        deriv[i+1] = (y[i+1]-y[i-1])/(x[i+1]-x[i-1])
    return deriv
    
def get_neutral_point(cg,Cm,CL):
    dCmCL = get_deriv(CL,Cm)
    NP = cg[0] - dCmCL
    return NP

def normalize_mac(x):
    xwing = 1.45
    mac = 1.2
    return (x-xwing)/mac

def get_thrust_mom_coef(cg,thrust):
    Sref = 11.4
    cref = 1.2
    V = 34.0
    rho = 1.2255
    #thrust = 130.0 #N
    dx = cg[2]-0.0
    return 2.0*thrust*dx/(rho*V*V*Sref*cref)

def output_report(alpha):
    pass

class OutResults:
    def __init__(self,ncase,alpha,name):
        self.name = name
        self.alpha = alpha
        self.n = len(alpha)
        self.ncase = ncase
        self.Cm = np.zeros([ncase,self.n])
        self.Cmref = np.zeros(self.n)
        self.SMref = np.zeros(self.n)
        self.SM = np.zeros([ncase,self.n])
        self.i = 0
    def add_rslt(self,CM,SM):
        self.Cm[self.i] = CM
        self.SM[self.i] = SM
        self.i += 1
    def show_rslt(self):
        h1, h2 = '',''
        for i in range(self.ncase):
            h1 += 'Cm %d\t'%(i+1)
            h2 += 'SM %d\t'%(i+1)
        header = 'alpha\tCm ref\t'+h1+'SM ref\t'+h2
        rslt = np.vstack([self.alpha,self.Cmref,self.Cm,self.SMref,self.SM])
        rslt = np.transpose(rslt)
        
        print '\n%s'%self.name
        print header
        for line in rslt:
            out =''
            for val in line:
                out += '%.4f\t'%val
            print out

def run_analysis():
    cgOld = [1.711,0.0,-0.104] #TODO: check values
    cgList = list()
    cgList.append([1.711,0.0,-0.104])
    cgList.append([1.711,0.0,-0.104])
    cgList.append([1.824,0.0,-0.125])
    cgList.append([1.824,0.0,-0.125])
    thrust = [1270.,0.,1270.0,0.]
    
    color = ['rs-','r--','bs-','b--']
    cgLegend = ['FWD power on','FWD power off','AFT power on','AFT power off']
    ncase = len(cgList)
    refLen = 1.2
    caseList = [9,10,14,15,18,19,22,23]
    caseTitle = ['AOS00,FLP00','AOS10,FLP00','AOS00,FLP10','AOS10,FLP10',
                 'AOS00,FLP25','AOS10,FLP25','AOS00,FLP35','AOS10,FLP35']
    
    for i in range(len(caseTitle)):
        caseTitle[i] = 'Case %d: '%(i+1) + caseTitle[i]
    
    for i,case in enumerate(caseList):
        wtRslt = wt.WTresult(case)
        output = OutResults(ncase,wtRslt.alpha,caseTitle[i])
        output.Cmref = wtRslt.CM
        fig = plt.figure()
        ax1 = fig.add_subplot(211)
        ax2 = fig.add_subplot(212)
        fig.suptitle(caseTitle[i])
        ax1.hold(True)
        ax1.grid(True)
        ax1.set_ylabel('Pitch moment')
        ax2.set_ylabel('Static margin')
        ax2.set_xlabel('AOA, deg')
        ax1.set_ylim([-0.7,0.3])
        ax2.set_ylim([-0.3,0.3])
        ax1.set_xlim([wtRslt.alpha[0],wtRslt.alpha[-1]])
        ax2.set_xlim([wtRslt.alpha[0],wtRslt.alpha[-1]])
        ax2.hold(True)
        ax2.grid(True)
        #ax1.plot(wtRslt.alpha,wtRslt.CM,'bs-')
        NPref = get_neutral_point(cgOld,wtRslt.CM,wtRslt.CL)
        SMref = (NPref-cgOld[0])/refLen
        output.SMref = SMref
        #ax2.plot(wtRslt.alpha,SMref,'bs-')
        box1 = ax1.get_position()
        box2 = ax2.get_position()
        ax1.set_position([box1.x0, box1.y0,box1.width*0.75,box1.height])
        ax2.set_position([box2.x0, box2.y0,box2.width*0.75,box2.height])
        for T,cg,c in zip(thrust,cgList,color):
            Cmthrust = get_thrust_mom_coef(cg,T)
            Cm = mrc_translate(cgOld,cg,refLen,wtRslt.CM,wtRslt.CL,wtRslt.CD,wtRslt.alpha) + Cmthrust
            NP = get_neutral_point(cg,Cm,wtRslt.CL)
            SM = (NP-cg[0])/refLen
            ax1.plot(wtRslt.alpha,Cm,c)
            ax2.plot(wtRslt.alpha,SM,c)
            output.add_rslt(Cm,SM)
        ax1.legend(cgLegend,loc='center left',bbox_to_anchor=(1,0.5))
        output.show_rslt()
    plt.show()


def debug1():
    refOld = np.array([1.711,0.0,-0.125])
    refNew = np.array([1.911,-0.029,-0.125])
    refLen = 1.2
    wtRslt = wt.WTresult(9)
    
    Cm = mrc_translate(refOld,refNew,refLen,wtRslt.CM,wtRslt.CL,wtRslt.CD,wtRslt.alpha)
    NP = get_neutral_point(refNew,Cm,wtRslt.CL)
    SM = (NP-refNew[0])/refLen

    plt.figure(1)
    plt.grid(True)
    plt.hold(True)
    plt.plot(wtRslt.alpha,Cm,'ro-')
    plt.plot(wtRslt.alpha,wtRslt.CM,'bo-')
    plt.show()


if __name__=="__main__":
    run_analysis()